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Delivery of anti-cancer drugs to tumor tissues, including their interstitial transport and cellu- 
lar uptake, is a complex process involving various biochemical, mechanical, and biophysical 
factors. Mathematical modeling provides a means through which to understand this com- 
plexity better, as well as to examine interactions between contributing components in a 
systematic way via computational simulations and quantitative analyses. In this review, we 
present the current state of mathematical modeling approaches that address phenomena 
related to drug delivery. We describe how various types of models were used to predict 
spatio-temporal distributions of drugs within the tumor tissue, to simulate different ways to 
overcome barriers to drug transport, or to optimize treatment schedules. Finally, we discuss 
how integration of mathematical modeling with experimental or clinical data can provide 
better tools to understand the drug delivery process, in particular to examine the specific 
tissue- or compound-related factors that limit drug penetration through tumors. Such tools 
will be important in designing new chemotherapy targets and optimal treatment strategies, 
as well as in developing non-invasive diagnosis to monitor treatment response and detect 
tumor recurrence. 
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INTRODUCTION 

Systemic chemotherapy is one of the most widely used treatments 
in all kinds of cancers and at every stage of tumor progression. 
However, success of the systemic treatment depends not only on 
the efficacy of chemical compounds, but also on whether these 
compounds can reach all tumor cells in concentrations sufficient 
to exert therapeutic effect. Most clinically used anti-cancer drugs, 
however, lead to the emergence of anti-drug resistance, and to over- 
come this therapeutic limitation, the chemotherapeutic agents are 
often used in combination with other drugs of different phar- 
macokinetic properties or in combination with other anti-cancer 
treatments. 

The process of drug delivery is complex and embraces dif- 
ferent temporal and spatial scales, including the organism level 
(where drug absorption, distribution, metabolism, excretion, and 
toxicity are studied in various organs and are known together 
under the acronym ADME-T), tissue and cell scales (where the 
main processes include drug extravasation into the tumor tis- 
sue, its penetration via interstitial transport, and cellular uptake), 
and intracellular level (where drug internalization, intracellular 
pharmacokinetics, accumulation, and efflux are investigated). In 
this review, we will focus on these mathematical models that act 
on the tissue scale. We refer the reader to the following research 
papers and review articles that address the other modeling scales 
(1-11). 



Transport of drug particles at the tissue level encounters several 
physiological and physical barriers. The architecture of tumor vas- 
culature is leaky and tortuous when compared to the vasculature 
of normal tissues. As a result, the blood flow is chaotic and the 
supply of nutrients and drugs irregular. This, in turn, leads to the 
emergence of regions of transient or permanent hypoxia. The cel- 
lular and stromal architecture of tumor tissue is far from being 
as well organized as that of normal tissues, and it is characterized 
by increased cell packing density, high variability in tumor cell 
sizes, and their locations. Together, these result in a non-uniform 
exposure of tumor cells to metabolites and drugs. Elevated inter- 
stitial fluid pressure (IFP), which is a consequence of the lack 
of functional lymphatic vessels, and vascular hyperpermeability, 
reduce extravasation of both fluid, and drug molecules from the 
vascular system, hindering advective transport through the tumor 
tissue. A dense extracellular matrix (ECM) with irregular align- 
ment of ECM fibers and with increased fiber cross-linking, also 
hinders the diffusion process. In general, it is difficult to pre- 
dict the extent of drug penetration into the tumor tissue and to 
determine the influence of various microenvironmental factors on 
drug interstitial transport. The former issue can be addressed by 
developing imaging techniques to visualize either the drug uptake 
or its lethal effects. The latter can be tested using systematical 
computational simulations of properly formulated mathematical 
models. 
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Several imaging approaches have been used to visualize the 
effects of drug penetration into the tumor tissue, including nat- 
urally fluorescent drugs showing their spatial distribution (12- 
14), specific imaging biomarkers showing the effects of anti- 
cancer drugs, such as cell DNA damage (15, 16), intravital 
microscopic imaging for real-time in vivo drug distribution (17), 
or molecular photoacoustic tomography (18). Numerous imag- 
ing techniques and their use in oncology have been reviewed 
in Weissleder and Pittet (19), Gillies et al. (20), and Morse and 
Gillies (21). 

Mathematical modeling provides tools for examining which of 
the various biophysical features of the tumor tissue and/or stroma 
and biochemical properties of drug compounds contribute sig- 
nificantly to limited drug penetration. In silico simulations are 
well-suited for testing combinations of multiple parameters that 
can be varied simultaneously in a controlled manner and over 
a wide range of values. Such a broad screening of drug or tis- 
sue conditions is rarely possible in laboratory experiments, but 
it is relatively easy and cheap in computer simulations. These 
theoretical screenings can help to determine the properties of ther- 
apeutic compounds optimal for their efficient interstitial transport 
(designing in silico drugs) or make decisions regarding the most 
effective drug combinations and scheduling protocols (design- 
ing in silico trials). Moreover, mathematical modeling allows for 
bridging laboratory experiments with clinical applications by pro- 
viding the means to extrapolate the in vivo results from mouse 
models to humans. Recently, several review papers discussing 
the power of mathematical and biophysical modeling have been 
published (22-29). 

In this review, we will focus on the most recent research articles 
that use mathematical and computational models of anti-cancer 
drugs acting on the cell and tissue scales. In the most general 
description, changes in the amount of drug present in the tis- 
sue depend on three values: the amount of drug entering the 
tissue (drug production), how the drug moves within the tissue 
(drug transport), and the amount leaving the tissue (drug elimi- 
nation). However, various phenomena can contribute to each of 
these three processes. For example, a drug can be supplied from 
the preexisting vascular system or can be released within the tis- 
sue from a moving drug carrier (such as a nanoparticle), or it 
can be activated due to specific environmental conditions (for 
example, low oxygen level or high acidity). Drugs can be carried 
through the tissue with the interstitial fluid flow (advective trans- 
port) or move randomly due to the Brownian motion of drug 
molecules (diffusive transport). Drug elimination from the tis- 
sue can take place due to its natural half-life (decay), binding to 
the ECM (degradation or deactivation), or cellular uptake. Math- 
ematically the simplest equation describing the kinetics of drug 
concentration c{x,t) at location x and at time t maybe written as 
follows: 



Here, k is a constant rate of drug supply, release, or activa- 
tion that takes place in a part of the domain (region), which may 
be a blood vessel (supply), nanoparticle (release), or low oxygen 
area (activation); D is a constant diffusion coefficient; u(x,t) is 
the velocity of the interstitial fluid; a is a decay or deactivation 
rate constant; and p is a rate constant of drug uptake by the cell. 
Schematically, all processes involved in the drug kinetics are shown 
in Figure 1 . Notably, each of these factors may take a more complex 
form. A more detailed discussion regarding these processes follows 
below, and we give examples of how they have been addressed in 
the mathematical modeling literature and applied to anti-cancer 
drug kinetics. 

MODELS ADDRESSING DRUG VASCULAR SUPPLY 

After intravenous infusion, drug molecules circulate in the vas- 
cular system before they extravasate into the surrounding tissue. 
The drug influx rate k is assumed constant in the equation listed 
above; however, more complex cases can be modeled wherein 
the vascular supply process depends not only on the molecule's 
size, but also on the physical properties of the vasculature and 
the target tissue. In general, small drug or metabolite molecules 
can cross the vascular wall more easily than larger molecules 
can, and they can extravasate into both healthy and tumorous 
tissues. Larger molecules, such as nanoparticles, require vascu- 
lar fenestration with larger pores to be able to leave the blood 
circulation system. Additional factors, such as electrostatic inter- 
actions between the particles and the negatively charged pores of 
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Intersitial Fluid Pressure (dimensionless) 




FIGURE 2 | Examples of typical outcomes from various 
mathematical models of drug penetration through the tumor 

tissue. (A) A hybrid model of tumor mass and discrete vasculature (left, 
red-tumor tissue, brown-vasculature) was used to investigate fluid and 
drug extravasation from the vasculature (right, color corresponds to 
drug concentration) [fromWu et al. (37), Figure 18]; (B) Patterns of 
diffusion-(top) and advection-(bottom) dominated transport of drug 
molecules (red dots) through the interstitial space between the cells 
(white circles) [from Rejniak et al. (58), Figure 6]; (C) A gradient of the 
interstitial fluid pressure (top) and transvascular pressure differences 



(bottom) in the growing tumor mass [from Stylianopoulos et al. (30), 
Figure 2]; (D) A drug concentration gradient from the vessel (white) 
outward with colors representing high (yellow) and low (brown) levels 
of diffusive drug [from Thurber et al. (17), Figure 8]; (E) Spatial 
distributions of hypoxic (yellow), necrotic (black), and apoptotic (green) 
tumor cells within the tumor mass treated with angiogenesis inhibitors 
for 37 weeks; the treatment is supplied from vasculature (red) [from 
Gevertz (34), Fig. 3]; (F) Structural adaptation of vessel diameters 
(colors represent the volume of blood flow) inside the tumors [from 
Pries et al. (31 ), Figure 6]. All figures reprinted with permissions. 



the vessel wall, have been studied by Stylianopoulos et al. (30). 
The mathematical model predictions suggested that electrostatic 
repulsion has a minor effect on the transvascular transport of 
nanoparticles, but electrostatic attraction, caused even by small 
cationic charges, can lead to a significant increase in the trans- 
vascular flux of nanoparticles into the tumor interstitial space 
(Figure 2C). 

The blood microcirculation within solid tumors is dysfunc- 
tional due to highly irregular vasculature (Figure 2F) that hinders 
delivery of both nutrients and drugs (31, 32). To investigate the 
distribution processes of small molecule drugs to cancer cells, a 
computational model based on fluorescent images of tumor func- 
tional vasculature was designed by Thurber et al. (17). The model 
was calibrated with experimental data and used to predict tem- 
poral changes in drug distribution profile around vessels with 
intermittent blood flow for a typical drug administration sched- 
ule (Figure 2D). Vascular images were also used by Baish et al. 
(33) to design a mathematical model that analyses drug diffu- 
sion in irregularly shaped domains based on two simple measures 
of vascular geometry. These include the maximum distance in 
the tissue from the nearest blood vessel and a measure of the 
shape of the spaces between vessels. This model can also pre- 
dict how new therapeutic agents that inhibit or stimulate vascular 



growth alter the functional efficiency of the vasculature within the 
tumor tissue. Computational simulations of vasculature-targeting 
agents and their influence on tumor growth have been also per- 
formed by Gevertz (34, 35). These biophysical models (Figure 2E) 
were used to explore the therapeutic effectiveness of two drugs 
that target the tumor vasculature, angiogenesis inhibitors (such 
as avastin) and vascular disrupting agents (such as combretas- 
tatin). The simulation results suggested that vasculature-targeting 
agents, as currently administered, cannot lead to cancer eradica- 
tion, although a highly efficacious agent may lead to long-term 
cancer control. The models, however, identified a treatment regi- 
men that can successfully halt simulated tumor growth, even after 
the cessation of therapy. 

Another computational study has been performed to test the 
effects that different drugs exert on the same mass of tumor tis- 
sue. Sinek et al. (36) compared the effectiveness of doxorubicin 
and cisplatin in vascularized tumors taking into account vascular 
and morphological heterogeneity. The simulation results showed 
that lesion-scale drug and nutrient distribution may significantly 
impact therapeutic efficacy. It has been also shown how the ther- 
apeutic effectiveness of doxorubicin penetration depends upon 
other determinants affecting drug distribution, such as cellular 
efflux and density, offering some insight into the conditions under 
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which otherwise promising therapies may fail and, more impor- 
tantly, when they will succeed. These simulations indicated that 
macroscopic environmental conditions, notably drug and nutrient 
distributions, give rise to considerable variation in lesion response, 
hence clinical resistance. Moreover, the synergy or antagonism 
of combined therapeutic strategies depends heavily upon this 
environment. 

The elevated IFP and high hydraulic conductivity can act like 
microenvironmental barriers for transvascular transport to both 
anti-cancer drugs and nutrients, as have been investigated by Wu et 
al. (37). It has been shown computationally that small blood vessel 
resistance and collapse may contribute to lower transcapillary flux 
of oxygen. Moreover, the higher IFP distribution in the simulated 
tumors affected oxygen extravasation negatively, which, in turn, 
hindered tumor growth by decreasing the oxygen transfer to the 
tissue (Figure 2A). In another study Pozrikidis (38) has investi- 
gated the overall hydrodynamics of the leakage problem through 
a permeable capillary, taking into account hydraulic conductivity 
of arterial, venous, and extravasation flow rates. This showed that 
interstitium dilation promoted the rate of extravasation. 

MODELS ADDRESSING DRUG RELEASE AND ACTIVATION 

To increase efficacy of therapeutic compounds and increase the 
time of drug survival inside the tumor tissue beyond its half-life, 
various methods of drug release and activation have been pro- 
posed. In our simple equation listed above, the release/activation 
rate k is defined as a constant, and the release/activation region is 
hypothetical. However, more complex mechanisms can be incor- 
porated in the models. The release region can represent a nanopar- 
ticle and can be varied in both space and time according to the 
changes in carrier locations; the activation rate may depend on 
local drug concentration or distribution of metabolites and may 
take place in hypoxic/acidic tumor areas, respectively, or may be 
stimulated by external factors, such as temperature or magnetic 
fields. 

Nanoparticles have gained much interest as potential carri- 
ers of therapeutic agents due to their size, which enable them 
to extravasate in the leaky tumor vasculature preferentially, and 
due to their modular functionality, which allows for release of the 
drug by controlled diffusion from the core across the polymeric 
membrane to the matrix. A mathematical model taking into con- 
sideration avascular tumor growth followed by angiogenesis and 
nanoparticle-based drug delivery has been applied by van de Ven et 
al. (39) to design optimal therapeutic protocols. In particular, the 
effects of nanoparticles carrying doxorubicin were simulated for 
various parameter values to determine how much drug per particle 
and how many particles need to be released within the vasculature 
to achieve remission of the tumor. Moreover, it has been shown 
that cell death on a population level is non-linear with respect 
to the drug concentration. The same team has simulated vascu- 
lar accumulation of blood-borne nanoparticles to analyze how 
nanoparticle vascular affinity depends on its size and ligand den- 
sity, as well as vascular receptor expression (40). It has been shown 
that for high vascular affinities, nanoparticles tend to accumu- 
late mostly at the inlet tumor vessels, leaving the inner and outer 
vasculature depleted of nanoparticles. For low vascular affini- 
ties, nanoparticles distribute quite uniformly in the intratumoral 



vasculature, but they exhibit low accumulation doses. It has been 
shown that an optimal vascular affinity can be identified by provid- 
ing the proper balance between accumulation dose and uniform 
spatial distribution of the nanoparticles. This balance depends 
on the stage of tumor development (vascularity and endothelial 
receptor expression) and the nanoparticle properties (size, lig- 
and density, and ligand-receptor molecular affinity). The timing 
and the location of drug release from nanoparticles have been 
investigated by Kim et al. (41) in a combination of in vitro exper- 
iments and mathematical modeling. It has been shown that gold 
nanoparticles carrying either fluorescein or doxorubicin molecules 
move and localize differently in an in vitro three-dimensional (3D) 
model of tumor tissue, depending on whether the nanoparticles 
are positively or negatively charged. Fluorescence microscopy and 
mathematical modeling show that uptake, not diffusion, is the 
dominant mechanism in particle delivery. These results indicate 
that positive particles may be more effective for drug delivery 
because they are taken up to a greater extent by proliferating cells. 
Negative particles, which diffuse more quickly, may perform better 
when delivering drugs deep into tissues. 

Another drug carrier, engineered macrophages, that are capa- 
ble of delivering pro-drugs to hypoxic areas within the tumor 
have been modeled by Webb et al. (42) and Owen et al. (43). 
In the former paper, two modes of action in the multicellular 
spheroids were investigated: either the macrophages delivered an 
enzyme that activated an externally applied pro-drug (bystander 
model), or they delivered cytotoxic factors directly (local model). 
The bystander model was comparable to traditional chemother- 
apy, with poor targeting of tumor cells in the center of the spheroid 
that are assumed hypoxic; on the other hand, the local model was 
more selective for the hypoxic regions. This work suggested that 
effective targeting of hypoxic tumor cells may require the use of 
drugs with limited mobility or whose action does not depend on 
cell proliferation. The latter article addressed a case where ther- 
apeutic macrophages were preloaded with nanomagnets and a 
magnetic field was applied to the tumor site. Both the conventional 
chemotherapy and chemotherapy with macrophages delivering 
hypoxia-inducible drugs were compared, and model simulations 
predicted that combining conventional and macrophage-based 
therapies would be synergistic, producing greater antitumor effects 
than the additive effects of each form of therapy. The model 
also revealed that timing is crucial in this combined approach 
with efficacy being greatest when the macrophage-based, hypoxia- 
targeted therapy is administered shortly before or concurrently 
with chemotherapy. 

The effects of applying heat to tumors treated with cisplatin 
have been investigated by El-Kareh and Secomb (44). A theoreti- 
cal model for the intraperitoneal delivery of cisplatin and heat to 
tumor metastases in tissues adjacent to the peritoneal cavity has 
shown increased cell uptake of drug, increased cell kill at a given 
level of intracellular drug, and decreased microvascular density. 
The model suggested that the experimental finding of elevated 
intracellular platinum levels up to a distance of 5 mm when the 
drug is delivered by a heated infusion solution is due to pene- 
tration of heat, which causes increased cell uptake of the drug. 
The effects of hyperthermia on chemotherapy were also investi- 
gated by Gasselhuber et al. (45) by developing a spatio-temporal 
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model of the release of doxorubicin from low temperature sen- 
sitive liposomes. This model showed that this treatment com- 
bined with thermal ablation allowed for localized drug delivery 
with higher concentrations in the tumor tissue than conventional 
chemotherapy. 

MODELS ADDRESSING DRUG DIFFUSIVE TRANSPORT 

In the model equation listed above, we used a constant diffusion 
rate D that leads to homogeneous diffusive transport. However, the 
diffusion may depend on the structure and other physical proper- 
ties of the tissue in which this process occurs. One extension of the 
above equation has been widely used in modeling the spread of 
gliomas in the brain where the diffusion in the white matter and 
gray matter was characterized by different diffusion coefficients 
(46, 47). 

In the context of drug penetration into the tumor tissue, 
Venkatasubramanian et al. (48) have created a mathematical model 
integrating intracellular metabolism, nutrient and drug diffusion, 
cell-cycle progression, and drug pharmacokinetics. Results indi- 
cated the existence of an optimum drug diffusion coefficient. A 
low diffusivity prevents effective penetration before the drug is 
cleared from the blood, and a high diffusivity limits drug reten- 
tion. This result suggests that increasing the molecular weight of 
the anti-cancer drug by nanoparticle conjugation would improve 
its efficacy. The simulations also showed that tumors that grow 
fast are less responsive to therapy than are tumors growing more 
slowly with greater numbers of quiescent cells, demonstrating the 
competing effects of regrowth and cytotoxicity. 

The complex interactions of drug particles and the ECM fibers 
that may hinder the drug molecule diffusion process have been 
modeled by Stylianopoulos et al. (49, 50). In this 3D model, sto- 
chastic fiber networks with varying degrees of alignment were con- 
sidered. Quantitative analysis of four different structures, ranging 
from nearly isotropic to perfectly aligned, were performed. The 
results indicated that the overall diffusion coefficient is not affected 
by the orientation of the network. However, structural anisotropy 
results in diffusion anisotropy, which becomes more significant 
with an increase in the degree of alignment, the size of the diffusing 
particles, and the fiber volume fraction. These model predic- 
tions were validated experimentally, showing for the first time 
in tumors that the structure and orientation of collagen fibers in 
the extracellular space leads to diffusion anisotropy. The authors 
also investigated the effects of charge on the diffusive transport 
of macromolecules and nanoparticles in the ECM, taking into 
account steric, hydrodynamic, and electrostatic interactions. The 
model showed that electrostatic forces between the fibers and 
the particles result in slowed diffusion. However, the repulsive 
forces become less important as the fiber diameter increases. These 
results suggest that optimal particles for delivery to tumors should 
be initially cationic to target the tumor vessels and then change to 
neutral charge after exiting the blood vessels. 

Since the ECM is composed of multiple cross-linked fibers, the 
drug particle diffusion in the interstitial space may rather resemble 
random movement through small nanochannels than diffusion 
through the open homogeneous space. A computational model 
that accounts for interface effects on diffusivity has been devel- 
oped and validated by predicting experimental glucose diffusion 



through a nanofluidic membrane (51-53). Moreover, the passive 
transport of nanoparticles from bulk into a nanochannel has been 
modeled, showing that subtle changes in nanochannel dimensions 
may alter the energy barrier. This results in different nanoparticle 
penetration depths and diffusion mechanisms. 

More detailed models of ECM structure, including fiber ori- 
entation, cross-link, and remodeling by the embedded cells have 
been developed by Bauer et al. (54, 55), Dallon and Sherratt (56) 
and Dallon et al. (57) in the context of vessel sprout and wound 
healing, respectively. These models have not yet been applied to 
model the role of ECM structure on drug molecule penetration. 
However, cellular heterogeneity of the stroma and its influence on 
both diffusive and advective forms of transport have been mod- 
eled by our group using idealized tissue morphologies of various 
porosity and cellularity values (58). Our simulations revealed that 
irregularities in the cell spatial configurations can solely result in 
the formation of interstitial corridors that are followed by drug 
or imaging agent molecules, leading to the emergence of tissue 
zones with less exposure to the drugs. Moreover, we showed that 
the relation between tissue porosity (defined as the extent of void 
space in the tissue), cellular density (defined as the number of cells 
per tissue area), and permeability (defined as time needed for a 
certain number of particles to traverse a predefined distance) is 
non-linear; thus it is also non-intuitive. 

MODELS ADDRESSING DRUG ADVECTIVE TRANSPORT 
During advective transport, drug molecules are carried with the 
flow of the interstitial fluid. This flow can arise from pressure dif- 
ferences within the tissue or from drainage of the fluid into the 
lymphatic circulation system. Wu et al. (37) investigated the role of 
the IFP, interstitial fluid flow, and the lymphatic drainage system on 
the transport of metabolites in developing tumors. The model sim- 
ulations showed that elevated interstitial hydraulic conductivity 
combined with poor lymphatic function is the root cause of the 
development of plateau profiles of the IFP in the tumor, which 
have been observed in experiments. 

At the macroscopic scale, where the individual cells are mod- 
eled as surrounded by the ECM space that is interpenetrated by the 
interstitial fluid, our group investigated the role of both advection 
and diffusion of drug molecules movement through the stroma 
(58). Simulation results collected from more than 100 different tis- 
sue morphologies showed that tissue cellular porosity and density 
influence the depth of drug penetration in a non-linear fash- 
ion. It has also been shown that for small diffusion coefficients, 
drug transport is advection dominated independently of tissue 
structure. Similarly, for all tissue structures considered in our sim- 
ulations, drug molecule transport was diffusion dominated for 
large diffusion coefficients. However, for the intermediate val- 
ues of fluid flow velocity and diffusion coefficients, the nature 
of interstitial transport depends strongly on the tissue morphol- 
ogy (Figure 2B). This indicates that sole knowledge of drug and 
tumor biophysical properties without knowledge of tumor tissue 
histology may lead to false predictions regarding the extent of drug 
penetration into the tumor tissue. 

The significant role of the advective fluid flow in brain tumors 
has been investigated by Arifin et al. (59, 60). In this work, a com- 
putational model was employed to simulate 3D patient-specific 
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distribution of carmustine. This model showed that a quasi-steady 
transport process is established within 1 day following treatment, 
and the drug is eliminated rapidly by transcapillary exchange, 
while its penetration into the tumor is mainly by diffusion. Con- 
vection appears to be crucial in influencing the drug distribution 
in the tumor resulting in non-homogeneous exposure to the drug: 
the remnant tumor near the ventricle is, by one to two orders of 
magnitude, less exposed to the drug than is the distal remnant 
tumor. In addition, local convective flow within the cavity appears 
to be a crucial factor in distributing the drug so that the tumor 
domain near the ventricle is prone to minimal drug exposure. The 
authors also simulated four chemotherapeutic agents (carmustine, 
paclitaxel, fluorouracil, and methotrexate) in a realistic 3D tissue 
geometry extracted from magnetic resonance images of a brain 
tumor. The simulation analysis showed that only paclitaxel exhib- 
ited minimal degradation within the cavity, as well as the best 
penetration of the remnant tumor. 

A mixture of computational modeling and laboratory exper- 
iments on gels and tumors reported in Ramanujan et al. (61) 
showed that the diffusive transport of drug particles might be 
obstructed more significantly by collagen fiber alignment than 
particle movement due to fluid advection. 

MODELS ADDRESSING DRUG DECAY. DEACTIVATION, AND 
CELLULAR UPTAKE 

In our simple equation above, the rate of drug decay, deactivation, 
and cellular uptake were defined as proportional to the local drug 
concentration. In the case of drug decay this is a typical way of 
incorporating drug half-life. In the case of drug deactivation or 
degradation, these processes may also depend on environmental 
factors, such as binding to ECM fibers or interacting with other 
mocroenvironmental factors. This aspect of drug pharmacody- 
namics has usually been neglected in mathematical models due to 
insufficient experimental data to inform or validate the models. 
However, with the recent advances in visualizing and experimen- 
tally quantifying ECM fibril structure, this process should be easier 
to incorporate in future mathematical models. Additionally, the 
process of cellular uptake can depend on various factors. Cer- 
tain drug molecules may bind to specific cell membrane receptors, 
and the efficacy of this process will then depend on the number 
of available receptors. Others may diffuse through the cell mem- 
brane, and this diffusion process will depend on both extracellular 
and intracellular drug concentrations. 

The complex interplay between molecular size, affinity, and 
tumor uptake has been investigated by Schmidt and Wittrup (62) 
using a mechanistic model that takes into account drug molecu- 
lar radius, interstitial diffusivity, available volume fraction, and 
plasma clearance. This model allowed for predicting the mag- 
nitude, specificity, time dependence, and affinity dependence of 
tumor uptake across a broad size spectrum of therapeutic agents. 
The authors concluded that the intermediate-size targeting agents 
(~25 kDa) have the lowest levels of tumor uptake, when compared 
to tumor uptake levels achieved by smaller and larger agents. In 
Thurber and Wittrup (63), this model was extended to create a 
mechanistic description of total antibody uptake in a tumor, tak- 
ing into account both free (unbound) antibody in the interstitium 
and antibody bound to its target. This allowed for an estimation 



of the time course of antibody uptake in solid tumors and its 
clearance from the blood plasma. 

The cellular pharmacodynamics of various anti-cancer drugs 
was investigated by a mathematical model that takes into account 
cellular uptake of the drug and both intracellular and extracel- 
lular cytotoxicities. In El-Kareh and Secomb (64), the damage 
induced by doxorubicin was expressed as the sum of two terms, 
representing the peak values over time of intracellular and extra- 
cellular drug concentrations. Drug uptake by cells was assumed to 
include both saturable and unsaturable components, which pro- 
vided better fits to in vitro cytotoxicity data. Model simulations 
suggested also a mechanism for the emergence of plateaus in the 
dose-response curve at high concentrations and short exposure 
time, as observed experimentally in some cases. Similar models 
were used to investigate the pharmacodynamics of cisplatin (65) 
and paclitaxel (66). 

TOWARD CLINICAL APPLICATIONS OF MATHEMATICAL 
MODELS 

Mathematical models can also provide the means to scale experi- 
mental results from animal to human body size and metabolism, 
and can be used to test various drug administration procedures 
and schedules (bolus injections, dose-dense therapies, continuous 
infusions, and adaptive therapies) in virtual human body. El-Kareh 
and Secomb (67) used mathematical modeling to determine the 
optimal mode of delivery for doxorubicin by comparing three 
intravenous administration methods: bolus injection, continuous 
infusion, and liposomal delivery. The model took into account the 
relatively slow rate and saturability of doxorubicin uptake by cells 
and predicted peak concentrations of drug attained in tumor cells, 
as well as peak concentration of free doxorubicin in blood plasma. 
The model simulations suggested that continuous infusion for 
optimal durations is superior to the other delivery methods. A 
similar model, but using the tumor cord geometry, was used by 
Eikenberry (68) to test doxorubicin dose optimization. Model sim- 
ulations showed that extending drug infusion time up to 2 h and 
fractionating large doses are two strategies that may preserve or 
increase anti-tumor activity, as well as reduce cardiotoxicity, by 
decreasing peak plasma concentration. Traina et al. (69) used the 
Norton-Simon tumor volume growth kinetic model (70) to pre- 
dict a tolerable dose of capecitabine (7 days treatment followed by a 
7-day rest) for advanced-stage breast cancer patients and this pre- 
diction was confirmed in phase I study. Traina et al. (7 1 ) continued 
to use the Norton-Simon model to optimize chemotherapeutic 
dosages and schedules in mouse xenograft models. Similar math- 
ematical models have been used to study dose-dense chemothera- 
pies (72) and to evaluate both the limitations of current schedules 
in breast cancer treatment and therapeutic advantages of novel 
dose-dense chemotherapies (73). Gatenby et al. (74) examined 
a novel approach in which cancer therapy was adapted to the 
evolving temporal and spatial variability of the tumor microenvi- 
ronment, cellular phenotypes, and therapy- induced perturbations 
instead of using a typical linear protocol of drug administration. 
The developed mathematical model suggested that if resistant pop- 
ulations are present before administration of therapy, the total 
elimination of the drug-sensitive subpopulation will lead to the 
faster growth of a drug-resistant population. As an alternative, 
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the simulated treatment was continuously modulated to control 
the size of tumor cell population that resulted in prolonged sur- 
vival. The authors went a step further and, actually, tested their 
predictions experimentally. In subsequent work, Silva et al. (3) 
parameterized the adaptive therapy model using in vitro exper- 
iments and showed that this treatment strategy delays tumor 
burden and increases time to progression in tumor models. 

The computational models are also well-suited to simulate 
treatments based on patient-specific parameters and tissue charac- 
teristics leading to personalized medicine. Our group investigates 
interstitial transport of drug and imaging agents using digitized 
samples of patients' tumor histology (58). Frieboes et al. (75) 
implemented a mathematical model of tumor drug-response that 
integrates simulations with biological data and includes the exper- 
imentally observed resistant phenotypes of individual cells. This 
integrative method could be used to predict resistance based on 
specific tumor properties, potentially improving treatment out- 
come. Kim et al. (76) uses a combination of micro- and macro- 
scopic imaging data and computational modeling to investigate 
blood flow in the heterogeneous tumor tissues. Venkatasubraman- 
ian et al. (77) uses breast cancer patients' DCE-MRI (dynamic 
contrast-enhanced magnetic resonance imaging) data to predict 
their responsiveness to therapeutic treatment. Their model sim- 
ulations showed that transvascular transport was correlated with 
tumor aggressiveness because of the formation of new vessels, 
and that increased transport heterogeneity led to increased tumor 
growth and poor drug-response. 

Clinically used imaging techniques will be crucial in integrating 
mathematical models with clinical data in order to make patient- 
specific predictions. For example, DCE-MRI technique allows for 
collecting time-activity curves with high spatial and temporal res- 
olution following a bolus injection of a Gadolinium-containing 
contrast agent, CA. The resulting data can be analyzed to generate 
spatially explicit (2D and 3D) maps of flow, perfusion, extra- 
cellular/ extravascular volume fraction and, in some cases, water 
permeability (78, 79). In general, the delivery and extravasation 
of CA is modeled with standard 2- or 3-compartment PK models 
(80) and, as mentioned above, these maps can be used to infer drug 
distribution in human tumors (77, 81). While many investigators 
use ROI (region of interest) analyses to derive a single perfusion 
value to describe a tumor, it is becoming increasingly appreciated 
that enhancement is heterogeneous and that quantitative descrip- 
tors of this heterogeneity improve the precision for diagnosis and 
monitoring of therapy response (82-86). We contend that per- 
fusion heterogeneity is a key factor in the response of tumors to 
therapy, both in terms of drug delivery and in the establishment of 
specific habitats that select for cells with specific phenotypes and 
hence, therapy responses (87). 

CONCLUSION AND FUTURE DIRECTIONS 

In this review, we discussed various mathematical models that 
were used to address different aspects of drug penetration through 
tumor tissue. All major stages of the penetration process have been 
investigated computationally: flow to different regions of tumors 
via blood vessels, crossing the vessel wall by drug molecules, 
their penetration through the interstitial tumor space, and cel- 
lular uptake. Mathematical models are well-suited to address such 



complex phenomena since by their nature they are able to handle 
multiple variables with numerous parameters. It is relatively easy 
and inexpensive to simulate tumor growth and treatment in sil- 
ico and to compare differences in simulation outcomes when such 
parameters are changed simultaneously and over a wide range of 
values. In fact, this area of mathematical research is dynamically 
expanding. Especially novel are models that account for spatial 
aspects of drug transport. Half of the papers described in this 
review and all of the images collected in Figure 2 come from man- 
uscripts that were published in the last 3 years, showing that this 
field is highly active and productive. 

Typically in mathematical models, the drugs are defined as con- 
centrations, as we did in the equation above. This is motivated 
by the fact that the number of drug molecules considered in the 
model might be a couple of orders of magnitude larger than the 
number of cells that form the tissue. However, in this description, 
only average behavior of drug molecules is captured. When more 
detailed drug kinetics need to be considered, such as molecule 
binding to cell receptors, intracellular trafficking, or mechanisms 
of drug extravasation from a vessel, drugs may be modeled as col- 
lections of individual molecules and can be traced individually in 
the model. Several novel models of this kind have been recently 
developed. Among models discussed in this review, the work of 
Ziemys et al. (51, 52), Mahadevan et al. (53), Frieboes et al. (40), 
and Rejniak et al. (58) traces the behavior of individual drug mol- 
ecules and their interactions with the cells and/or vessels. In the 
first two papers the authors also discuss how to scale between the 
description of the kinetics of individual drug molecules and more 
general description of drug concentration. 

Similarly, the more classical modeling approaches consider 
tumors as large populations of cells and represent them as cell den- 
sities (25, 36, 37, 40, 59, 64, 67, 68, 71, 73, 74, 81). These models can 
handle multiple cell subpopulations, but the number of different 
cell types has to be defined a priori, and new subpopulations can- 
not be dynamically created during the simulation. However, under 
specific conditions cells can be moved from one subpopulation to 
another. The predefined cell types may include a specific phase of 
the cell-cycle (a population of cells in Gl, GO, S, or G2 phase), 
a particular cell phenotype (a population of proliferating, quies- 
cent, hypoxic, or necrotic cells), or a particular cell response to 
the treatment (a population of drug-resistant or drug-sensitive 
cells). The advantage of continuous models is that they can han- 
dle large populations of cells, but the significant disadvantage is 
that all cell properties in these models must be averaged, since no 
individual cells are considered. In view of the growing evidence 
of heterogeneity of tumor cells on the genetic, phenotypic, and 
drug-response levels, the averaged cell properties and the averaged 
cell responses to anti-cancer treatments may not be sufficient to 
make predictions for individual patients. 

In contrast, in the individual-cell-based models (called also the 
single-cell-based models, or the agent-based models) each cell 
is represented as a separate entity that acts as an independent 
agent according to some predefined rules (cell phenotype), but 
cell behavior can also be modulated by interactions with other 
cells and with the immediate cell microenvironment (selection 
forces). In this class of models cells may differ from each other 
significantly (cells may have distinct phenotypes, independently 
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regulated cell-cycles, different levels of receptors, or different accu- 
mulations of mutations). Several models discussed in this review 
are single-cell-based (34, 35, 43, 52, 53, 58). The main advantage 
of these models is their natural cellular heterogeneity that better 
represents tumor multicellular composition than the continuous 
models. It is, of course, possible to analyze results of individual- 
cell-based models on a cell population level (in terms of average 
values and standard deviations, distributions, or correlations), 
similarly as this is done with experimental measurements. More- 
over, these analyses can be compared to results from continuous 
models. The inverse process, that is extracting detailed information 
on individual cells from continuous models, is impossible. The 
main disadvantage of agent-based models is in their limitations 
to handle large number of cells. Typically, this limit is in thou- 
sands of cells, but with constantly increasing speed of computers 
and with development of novel faster computational techniques 
(parallel, GPU, or cloud computing) the number of cells that the 
model can handle in a reasonable time may not be a limitation 
anymore. 

It is worth noting that every mathematical model is by its nature 
a simplification of the biological system it is assumed to represent, 
so we do not expect that one model will incorporate all processes 
involved in drug penetration through the tumor tissue. And we 
also do not expect that the unified modeling framework addressing 
all aspects of drug transport through tumor tissue will emerge in 
the near future. In silico models need to be designed to investigate 
a specific research question similarly to how biological experi- 
ments focus on the selected aspects of tumor treatment and do 
not address in a single experiment all possible combinations of 
involved factors. Computational models should not be too com- 
plex to allow for quantitative analysis of the relative importance 
of all features and parameters included in the model. However, in 
contrast to experiments, model parameters (e.g., drug molecular 
mass or charge, timing and dosing of drugs, and their activation 
or uptake properties) can be varied over a wide range of values 
and can be changed simultaneously in a controlled way, giving 
investigators insight into a full spectrum of drug properties that 
lead to the desired (or undesired) effects. These model outcomes 
will then provide guidance for further laboratory experimenta- 
tion, and both results, positive and negative, will be informative 
for biologists. The positive results will suggest the environmen- 
tal conditions or drug concentrations that are worth pursuing 
experimentally; the negative results will advise the drug concen- 
trations or their properties that do not lead to a desired effect 
and can be omitted, reducing experimental costs and time. In 
fact, close collaboration between mathematical modelers, biolo- 
gists, and clinicians is crucial, in our opinion, for making progress 
in improving anti-cancer treatments. 

In our opinion the computational models of tumor develop- 
ment and treatment that will be successfully applied in personal- 
ized medicine need to be single-cell-based to be able to account 
for differences between tumors in individual patients (inter-tumor 
heterogeneity) and between distinct regions within the same 
tumor tissue (intra-tumor heterogeneity) . Such models will be able 
to address phenotypic, genetic, and drug-response heterogeneity 
observable in patient tumors. The future models need to be tempo- 
ral to capture the dynamics of tumor growth, cell-cell interactions, 



and response to therapy. These models will allow for temporal 
analysis of model results in order to identify more effective drug 
administration schedules with potentially variable schedules and 
dosages that cannot be intuitively inferred from analyzing drug 
properties in laboratory experiments. The future models should 
also be spatially explicit and three-dimensional, since both cell 
growth dynamics and drug transport dynamics are significantly 
different between the one-, two-, and three-dimensional spaces. 
In vivo tumors have complex geometries, variable cellular densi- 
ties, irregular vasculature that cannot be captures by simple non- 
spatial, population-based models. And over all the future models 
need to be quantitative, based on quantitative experimental data 
(to inform and parameterize the model), and producing quantita- 
tive results, that can be compared to experimental measurements 
or clinical data. 

Given the complexity of processes taking place during tumor 
development and its treatment, as well as significant inter-patient 
and intra-tumor variability, the cross-disciplinary approaches that 
integrate data and methods from various scientific disciplines have 
a better chance to delineate the mechanisms of tumor resistance 
to treatment and the way to overcome drug delivery barriers. The 
mathematical models that are properly integrated with experimen- 
tal data, such that both in silico models and laboratory experiments 
inform each other, can provide tools for interpreting data, evaluat- 
ing the most important parameters for designing new experiments, 
and developing strategies to improve tumor treatment. 
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